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Abstract 



Computational geometry concepts are often easiest to understand visually, in terms 
of the geometric objects they manipulate. Indeed, most papers in the field rely on 
diagrams to communicate the intuition behind their results. However, static figures 
are not always adequate. 

The accompanying videotape showcases advances in the use of algorithm ani- 
mation, visualization, and interactive computing in the study of computational ge- 
ometry. This report contains brief descriptions of all the segments of the videotape. 

This report is the third annual Video Review of Computational Geometry to be 
published as a SRC Research Report. The '92 Video Review is SRC Research Re- 
port 87, and the '93 Video Review is SRC Research Report 101. Future editions 
will be published by the ACM, as part of the proceedings of the annual ACM Sym- 
posium on Computational Geometry. 



Preface 



Computational geometry concepts are often easiest to understand visually. Indeed, 
most papers in computational geometry rely on diagrams to communicate the in- 
tuition behind their results. However, static figures are not always adequate to de- 
scribe geometric algorithms, since algorithms are inherently dynamic. The accom- 
panying videotape showcases advances in the use of algorithm animation, visual- 
ization, and interactive computing in the study of computational geometry. 

This report contains short descriptions of the eight segments in the videotape. 
The first three segments are animations of two-dimensional algorithms: a polygonal 
approximation algorithm from cartography, an algorithm for computing the rectan- 
gle discrepancy, and a fixed-radius neighbor-searching algorithm. The fourth seg- 
ment presents the GASP system for animating geometric algorithms in two or three 
dimensions; GASP is used by two other segments in the review. The remaining 
segments illustrate algorithms and concepts in three or more dimensions: the fifth 
segment explores the relationship between penumbral shadows in three dimensions 
and four-dimensional poly topes; the final three segments show three-dimensional 
algorithms for interactive collision detection, polyhedral separators, and interactive 
display of alpha hulls. 

This report and the accompanying videotape comprise the Third Annual Video 
Review of Computational Geometry. In the three years since the Video Review was 
introduced, it has become an established part of the annual ACM Symposium on 
Computational Geometry. The 1995 Video Review and future editions will be dis- 
tributed by the ACM as part of the conference proceedings. 

We thank Bob Taylor, director of the Systems Research Center, for his support 
in establishing the annual Video Review. The '92 Video Review is available as SRC 
Research Report 87, and the '93 Video Review as SRC Research Report 101. These 
reports (written guides and videotapes) are available by sending electronic mail to 
src-report@pa . dec . com. 

Finally, we thank the members of the 1994 Video Program Committee for their 
help in evaluating the entries. The members of the committee were Marshall Bern 
(Xerox PARC), Marc Brown (Digital Equipment Corporation's Systems Research 
Center), Leo Guibas (Stanford University), John Hershberger (Mentor Graphics), 
Jim Ruppert (NASA Ames), and Jenny Zhao (Silicon Graphics). 

Marc H. Brown 
John Hershberger 
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The "line simplification problem" is as follows: Given a sequence of vertices that 
make up a polygonal chain, find a subsequence of the vertices that approximates 
the chain. The accompanying videotape shows two algorithms for solving the line 
simplification problem. 

Introduction 

An important task of the cartographer's art is to extract features from detailed data 
and represent them on a simple and readable map. As computers become increas- 
ingly involved in automated cartography, efficient algorithms are needed for the 
tasks of extraction and simplification. Cartographers [1,6] have identified the line 
simplification problem as an important part of representing linear features. Given 
a polygonal chain, a sequence of vertices V = {i>o, v\, . . . ,v n }, the problem is to 
find a subsequence of vertices that approximates the chain V. We assume that the 
input chain V has no self-intersections (but not the output). 

In 1973, Douglas and Peucker (now Poiker) [3] published a simple recursive 
simplification method in the cartography literature that has become the favorite of 
many [5,10]. This method has been independently proposed in other contexts such 
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Figure 1: Simplifying the polygonal chain representing the 
coastline of British Columbia. 



as vision [8] and computational geometry [9]. The accompanying videotape ani- 
mates two different implementations of this simplification method — the straightfor- 
ward implementation suggested by Douglas and Peucker and the implementation of 
Hershberger and Snoeyink that improves the worst-case running time. 

Method 

Given a polygonal chain V = {i>o, v\, . . . ,v n }, the recursive line simplification 
method proceeds as follows: Initially, approximate V by the line segment vov n . 
Determine the farthest vertex vj from the line 1)qv^. If its distance d(vj, IJqv^,) is 
at most a given tolerance e > 0, then accept the segment vqv^ as a good approxi- 
mation to V. Otherwise, break V at Vf and recursively approximate the subchain 
{d 0 , vi,...,Vf} and {v f , v n }. 

Algorithms 

The two implementations differ in the way they find Vf, the farthest vertex from 
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The straightforward approach measures distance from^o^ for each vertex and 
takes the maximum. Thus, the running time of this algorithm will satisfy a quicksort- 
like recurrence, with best case G(n log n). Since geometry determines the farthest 
vertex, the worst-case running time is 0(n 2 ). One cannot use linear- time median 
finding or randomization to improve the running time. (If one assumes that the in- 
put data is such that farthest vertices are found near the middle of chains, then one 
expects O(nlogn) behavior.) 

The second implementation uses a path hull data structure [2, 4] to maintain a 
dynamic convex hull of the polygonal chain V. Using Melkman's convex hull al- 
gorithm [7], we compute two convex hulls from the middle of the chain outward. 
We can find a farthest vertex v / from a line by locating two extreme points on each 
hull using binary search. When we split V at vertex vj, we undo one of the hull 
computations to obtain the hull from the "middle" to vj, recursively approximate 
the subchain containing the middle, then build hulls for the other subchain and re- 
cursively approximate it. One can use a simple credit argument to show that the 
total time taken by the algorithm is 0(n log n). The best case can even be linear, if 
all the hulls are small and splits occur at the ends. 

Unfortunately for us, the statistical properties of cartographic data usually mean 
that the straightforward implementation is slightly faster than the path hull imple- 
mentation. Since the variance of the running time of the path hull implementation 
is smaller, it may be interesting for parallel or interactive applications. 

Technical Details 

The programs were run on a Silicon Graphics Crimson in the GraFiC lab at UBC 
and output was recorded to video in real time. 

Acknowledgments 
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Computing the Rectangle Discrepancy 

David P. Dobkin and Dimitrios Gunopulos 
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Supersampling is one of the most general ways of attacking the problem of antialias- 
ing in computer graphics. In this approach, a picture is sampled at a very high rate, 
and then resampled by averaging supersamples within the pixel area to compute the 
pixel value. Ray tracing, a common technique used to produce realistic images of 
computer modeled scenes, is another instance of supersampling. Discrepancy the- 
ory provides a measure of the quality of the sampling patterns. The "rectangle dis- 
crepancy" in particular gives a good measure of how well a sample pattern, when 
applied to the area of one or more pixels, captures small details in the picture. The 
accompanying videotape shows an animation of an algorithm we have developed 
[2] for computing the "rectangle discrepancy" of a two-dimensional point set. In 
particular, this algorithm can be used to find sampling sets with very low rectangle 
discrepancy. Interestingly the problem of computing the rectangle discrepancy also 
arises in the context of learning theory [4]. 

Definitions 

Let us consider a point set X in the unit square. Let T be the set of all rectangles in 
the unit square with edges parallel to the two axes. For an axis-oriented rectangle 
A (A € T), we define the function 

In{A) = \AC\X\ 
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max (In(A)/IXI - Area(A)) 




Figure 2: The rectangle that maximizes T>i. 

The rectangle discrepancy of a rectangle A G T is 

V{A) = \In(A)/\X\ - Area(A)\ 
The maximum rectangle discrepancy of the set X is 

MaxV(X) = max Ae r(V(A)) 
We also define the following simpler functions for A £ J 7 

V l (A) = In{A)/\X\ - Area{A) 

and 

V 0 {A) = Area{A) - In{A)/\X\ 
The algorithm we present in this animation computes 

MaxVi(X) = max Ae r(Vi(A)) 

for a given input set X, in time O ( | X | 2 log | X | ) . The same techniques however can 
be used to maximize V 0 , and so this problem is equivalent to finding MaxV(X). 



6 



2 



Videotape Contents 

We explain the main idea of the algorithm in a series of simple lemmata. First we 
show that we have to consider at most Od^ 4 ) rectangles, more specifically the 
ones that have all their edges pass through a point in X. 

This allows us to study pairs of points in in X. For each such pair we find 
the rectangle that maximizes V\ over all rectangles with horizontal edges that pass 
through the points of the given pair. 

This problem is much easier now because we know the y coordinates of the 
rectangle. Assuming that the two points are pb,Pt (Pb{y) < Pt{y)X and the height 
of the rectangle is dy, (dy = pt(y) — Pb{y)) we have to find the following maximum 

max AeTij (In(A) j\X\ - dx A dy) 

where Tij is the set of rectangles that pass through points pb and p t , and dx A is the 
width of A. If we eliminate the points in X that either lie above the top y coordinate 
of the rectangle or below the lower y coordinate, this becomes a one-dimensional 
problem. We project the points on the x axis and obtain a new set X' of one-dimensional 
points. We can equivalently maximize the following function 

max Ae yQ^ Ph{x)eA {In{A)l\X'\ - Length{A)) 

We can reduce this problem to the problem of computing the halfspace discrepancy 
of one-dimensional point sets, a simpler problem studied in [1]. Here we want to 
find the anchored interval that maximizes 

max A=[0}Zi] (i/\Z\ - z t ) 

for a given one-dimensional set Z = {z\ . . . z n } with 0 < Zi < Zi + \ < 1. This 
is a linear function of the rank and the coordinate of the point. To solve it we find 
the convex hull of the set {(z\, 1) . . . (z n , n)} and find the maximum with a binary 
search. From this maximum we obtain the solution to the original maximization 
problem. 

In this way we compute the maximum rectangle for all pairs of points. Clearly 
the rectangle that maximizes V\ is the maximum of these 0(\ X\ 2 ) rectangles. 

Finally, we show how to use the dynamic data structure of [5] to solve a se- 
quence of one-dimensional problems efficiently. With the use of this technique we 
can maintain the convex hull at a cost of 0(log 2 \X\) time per insertion. This in 
turn allows us to find the maximum rectangle for each pair of points in O (log 2 | X \ ) 
time, and yields the final running time of 0(\X | 2 log 2 \X\). 
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Technical Details 

The accompanying videotape was prepared in the Computer Science Department 
at Princeton University. The implementation of the algorithm was written in C, to 
ran under UNIX on a Silicon Graphics Iris workstation. The animation was created 
with the animation system GASP that Ayellet Tal and David Dobkin are developing 
in this department [3]. Recording was done at the Interactive Computer Graphics 
Lab at Princeton and editing was done with the assistance of the Department of Me- 
dia Services at Princeton. 

Acknowledgments 
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An Animation of a Fixed-Radius All-Nearest 
Neighbors Algorithm 

Hans-Peter Lenhof and Michiel Smid 
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A fixed-radius all-nearest neighbors problem takes as input a set S of n points and 
a real number 6, and reports all pairs of points that are at distance at most 6. The 
accompanying videotape shows an animation of an algorithm, developed by the au- 
thors [2], for solving the fixed-radius all-nearest neighbors problem. 

The Algorithm 

The algorithm depicted in the videotape works in two stages. In the first stage, a 
grid is computed. Instead of a standard grid, we use a so-called degraded grid that 
is easier to construct by means of a simple sweep algorithm. This degraded grid 
consists of boxes with sides of length at least 6 . If a box contains points of S, then 
its sides are of length exactly 8. In the second stage, this grid is used for the actual 
enumeration. For each non-empty box, it suffices to compare each of its points with 
all points in the same box or in one of the neighboring boxes. Although the algo- 
rithm may compare many pairs having distance more than 6, it can be shown that 
the total number of pairs considered is proportional to the number of pairs that are 
at most 6 apart. As a result, the total running time of the algorithm is proportional 
to n log n plus the number of pairs that are at distance at most 6. 

As an application, we give an animation of the algorithm of Heiden et al. [1] for 
triangulating the contact surface of a molecule. In a first step, points on this surface 
are computed. Given these points, we compute all pairs that are at distance at most 
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8 for a suitable choice of 8. This gives for each point a neighbor list containing all 
points that are at distance at most 8 from it. Given these lists, we can triangulate 
the surface incrementally: An edge of the current triangulation is called an outer 
edge, if it belongs to only one triangle. All outer edges are maintained in a queue. 
With each edge, we store a so-called common neighbor list, which is the intersection 
of the two neighbor lists of the points belonging to this edge. We start with any 
triangle. Its edges and their common neighbor lists are inserted into the queue. As 
long as the queue is non-empty, we do the following: We take any edge from the 
queue and look at all points in its common neighbor list. Each such point defines a 
triangle with the current edge. We take that point defining the smallest triangle. This 
triangle becomes part of the triangulation and we update the queue appropriately. 

Technical Details 

The video was produced at the Max-Planck-Institute for Computer Science, using 
Silicon Graphics Indy and Indigo 2 machines, and using Showcase and Geomview- 
1.1 software. We thank Roland Berberich, Jack Hazboun, Dirk Louis and Peter 
Miiller for their help in producing this video. 

Acknowledgments 

This work was supported by the ESPRIT Basic Research Actions Program, under 
contract No. 7141 (project ALCOM II). 
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The accompanying videotape shows five algorithm animations that have been de- 
veloped using GASP. The GASP system is noteworthy because it is very easy for a 
user without any knowledge of computer graphics to create an animation of a geo- 
metric algorithm quickly. 



Overview 

The GASP system helps in the creation and viewing of animations for geometric 
algorithms. The user need not have any knowledge of computer graphics in order 
to quickly generate an animation. GASP allows the fast prototyping of algorithm 
animations. A typical animation can be produced in a matter of days or even hours. 
Even highly complex geometric algorithms can be animated with ease. The system 
is also intended to facilitate the task of implementing and debugging geometric al- 
gorithms. 

The application's code specifies only the structure of the animation (which build- 
ing blocks are included), and the structure of each scene (e.g., the position of the 
objects). The user need not be concerned with the way the animation appears on 
the screen. If the user wishes to influence the look of the animation and not only 
its structure, he can edit an ASCII "Style File" which controls viewing aspects of 
the animation. Even when the user changes the style file, the user needs no specific 
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Figure 3: Polyhedral Hierarchies 



knowledge of computer graphics. The animation is generated automatically by the 
system. But a different animation will be generated if the style file is modified. 

The GASP environment allows the viewer to control the execution of the al- 
gorithm in an easy way. The animation can be run at varying speeds: fast(>>), 
slow(>) or step by step (> |). The analogous <, << and | < push buttons run 
the algorithm in reverse. The viewer can PA USE at any time to suspend the execu- 
tion of the algorithm or can EJECT the movie. The environment is invaluable for 
education and debugging. 

Videotape Contents 

This video shows five examples of GASP in action. These examples were chosen 
to give a sense of the range of possible applications of GASP. In addition, two other 
videos [1], [3] that use our system appear in this collection of animations. 
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Clip 1: Polyhedral Hierarchies 

To begin, we show a clip from [6]. This is the hierarchical data structure for rep- 
resenting convex polyhedra [4, 5], and the animation was an aid in debugging the 
algorithm for detecting plane-polyhedral intersection. We explain a single pluck 
and then show how the hierarchy progresses from level to level. To create the an- 
imation for explaining a single pluck the user need only write the following short 
code: 

explain_pluck (poly_vert_no, poly_vertices , 
poly_f ace_no, poly_faces, poly_names , 
vert_no, vertices, face_no, faces) 

{ 

/* create the polyhedron */ 
Begin_atomic ( "poly" ) ; 
Create_polyhedron ("PO", 

poly_vert_no, poly_f ace_no, poly_vertices , poly_f aces ) ; 
Rotate_world ( ) ; 
End_atomic ( ) ; 

/* remove vertices and cones */ 
Begin_atomic ( "pluck" ) ; 

Split_polyhedron (poly_names , "PO", vert_no, vertices); 
End_atomic ( ) ; 

/* add new faces */ 

Begin_atomic ( "add_f aces" ) ; 

Add_f aces (poly_names [ 0 ] , face_no, faces); 

End_atomic ( ) ; 

/* undo plucking */ 
Undo ( 2 ) ; 

} 
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Clip 2: Objects That Cannot Be Taken Apart With Two Hands 

The second animation, based on [10], is a remake of [9]. The animation shows a 
configuration of six tetrahedra that cannot be taken apart by translation with two 
hands. Then, it presents a configuration of thirty objects that cannot be taken apart 
by applying an isometry to any proper subset. The purpose of the animation is to 
illustrate the use of GASP as an illustration tool for geometric configurations. It 
took us far less than a day to generate that animation. The animation was created 
by the following function: 

hands ( ) { 

for (i =0; i < stick_no; i++) { 
/* stick i */ 
get_polyhedron (spoints, 

Sindices, Snmax, Sfmax, &atomic_name, &stick_name) ; 

Begin_atomic (atomic_name) ; 

Create_polyhedron (stick_name, nmax, fmax, points, indices); 
End_atomic ( ) ; 

} 

Begin_atomic ( "Rotate " ) ; 
Rotate_world ( ) ; 
End_atomic ( ) ; 

} 

The following is part of the style file for the above animation. The style file 
determines the following aspects of the animation. The background color is light 
gray. The colors to be chosen by GASP are colors which fit the creation of a video 
(rather than the screen). Each atomic unit spans 10 frames. Rotation of the world is 
done 380 degrees around the Y axis, and it is executed over 160 frames, instead of 
over 10. The color of the stick which fades in during atomic unit stickO is red. 
(Similarly, the style file includes colors for the other sticks.) 
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begin_global_sty 
background = 0 
color = VIDEO; 
frames = 10; 

end_global_style 

begin_unit_style 
frames = 160; 
rotation_world 

end_unit_style 

begin_unit_style 
color = 0 . 625 ' 
end_unit_style 



9 0.9 0.9; 

Rotate 

= Y 380.0; 

stickO 

.273475 0.273475; 
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Clip 3: Motion Planning 

Next we show an animation of a motion planning algorithm [7]. It is a repeat of 
the first part of [8]. The object to be moved is a disc and the obstacles are simple 
disjoint polygons. The goal of the animation is to show the use of GASP in creating 
two-dimensional animations. 

Clip 4: Line Segment Intersection 

The fourth animation, which is based on [2], shows a line segment intersection al- 
gorithm in action and illustrates its most important features. The animation consists 
of three phases. The first phase presents the initial line segments and the visibility 
map that needs to be built. The second phase demonstrates that the visibility map is 
being constructed by operating in a sweepline fashion, scanning the segments from 
left to right, and maintaining the vertical visibility map of the region swept along 
the way. The third phase demonstrates that the cross section along the sweepline is 
maintained in a lazy fashion. The animation is used as an aid in explaining a highly 
complex algorithm and to allow students to interact with the algorithm. A student 
can choose the initial line segments by editing an ASCII file and view how the an- 
imation is changing. 

Clip 5: Heapsort 

The last clip animates heapsort. The colored rods have lengths to be sorted. We first 
display the initial creation of the heap. Next, simple motions are used to control 
the repeated operations of removing the largest element, and reheaping. Heapsort 
is given as an example of using GASP to facilitate the animation of any algorithm 
that involves the display of three dimensional geometries. 

Technical Details 

This video was prepared in the Computer Science Department at Princeton Uni- 
versity. The programs were written in C to run under UNIX on a Silicon Graphics 
Iris. Recording was done at the Interactive Computer Graphics Lab at Princeton 
and editing was done with the assistance of the Princeton Department of Media Ser- 
vices. 
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The shadows cast by polygonal occluding objects illuminated by polygonal light 
sources are complicated. Unlike shadows cast by point sources, they are divided 
into regions where the shadow intensity changes at different rates. The boundaries 
between these different regions form the shadow diagram. 

The accompanying videotape depicts shadows cast by a square light source with 
a square occluding object. We have included two types of images: flat shadow 
scenes and perspective views. The flat shadow images were calculated using the 
exact formula for radiative energy transfer. This calculation depends only on the 
boundary of the visible light source. For polygons, it can be reduced to a simple 
summation [1]. 

Finding the boundary of the visible portion of the illuminating polygon when it 
is occluded by a second polygon requires computing the intersection of two poly- 
gons. General purpose algorithms for computing intersections of arbitrary polygons 
in the plane are quite complicated, so we used the Sutherland-Hodgman algorithm, 
which cannot handle concave polygons, but which is much simpler than the fully 
general methods [2]. We created the three-dimensional perspective views from the 
flat images using the texture mapping capabilities of Rayshade. 

The shadow diagrams can be obtained by projecting the edges or vertices of the 
light source through vertices or edges of the occluder onto the plane of the shadow. 
The video illustrates that when a square occluding object is illuminated by a square 
light source, the shadow diagram is a projection of a four-dimensional hypercube. 
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In general, the shadow diagrams for an arbitrary light source and occluder can 
be seen as projections of four-dimensional polytopes by considering the lines which 
intersect the plane of the shadow. Let S be the plane of the shadow and let P be a 
different plane parallel to S. If a line intersects S at the point (x s , y s ) and it inter- 
sects P at the point (x p , y p ) then the line can be mapped to a point (x s , y s ,x p ,y p ) 
in R 4 . Now let D be the set of all lines which intersect both the light source and the 
occluder. When D is mapped into R 4 , the result is a four-dimensional poly tope. 
If this polytope is projected into R 2 by discarding the last two coordinates of each 
point, the shadow diagram results. Because P can be changed, this construction 
describes a particular shadow diagram as the projections of any one of a family of 
polytopes. 
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The accompanying videotape demonstrates algorithms for collision detection in large- 
scaled, interactive environments. Such environments are characterized by the num- 
ber of objects undergoing rigid motion and the complexity of the models. The algo- 
rithms do not assume that the motions of the objects are expressible as closed form 
functions of time, nor do they assume any upper bounds on their velocities. 

Background 

Over the last few years, a great deal of interest has been generated in virtual or 
synthetic environments, such as architectural walkthroughs, visual simulation sys- 
tems for designing the shapes of mechanical parts, training systems, etc. Interac- 
tive, large-scaled virtual environments pose new challenges to the collision detec- 
tion problem, because they require performance at interactive rates for thousands 
of pairwise intersection tests [1]. 

Our algorithms use a two-level hierarchical detection test for each model to se- 
lectively compute the precise contact between objects, achieving real-time perfor- 
mance without compromising accuracy At the top level, we use a dimension reduc- 
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Figure 5: An example of collision detection for convex poly- 
hedra 

tion approach (based on a technique termed as sweep and prune) and the temporal 
and geometric coherence that exists between successive frames to overcome the 
bottleneck of 0(n 2 ) pairwise comparison tests in a large environment of n moving 
objects [1]. After we eliminate the object pairs that are not in the close vicinity of 
each other, we employ an exact collision detection scheme. The exact convex poly- 
tope collision detection algorithm is based upon [2] : locality (captured by walking 
from one Voronoi region of the candidate feature to one of its neighboring features' 
Voronoi regions according to the constraints failed) and coherence (used to reduce 
computation efforts in a dynamic environment). We can deal with non-convex ob- 
jects or articulated bodies by using a sub-part hierarchical tree representation. 

Exact Collision Detection 

The Voronoi regions form a partition of space outside the polytope according to the 
closest feature. The collection of Voronoi regions of each polytope is the general- 
ized Voronoi diagram of the polytope. Note that the generalized Voronoi diagram of 
a convex polytope has linear size and consists of polyhedral regions. Each Voronoi 
region is bounded by a set of constraint planes with pointers to the neighboring cells 
(which each share a constraint plane with it). 

If a point lies on a constraint plane, then it is equidistant from the two features 
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which share this constraint plane in their Voronoi regions. If a point P on object 
A lies inside the Voronoi region of the feature on object B, then is a closest 
feature to the point P. 

Our method is straightforward. We start with a candidate pair of features, one 
from each polytope, and check whether the closest points lie on these features. Since 
the polytopes and their faces are convex, this is a local test involving only the neigh- 
boring features of the current candidate features. If either feature fails the test, we 
step to a neighboring feature of one or both candidates, and try again. With some 
simple preprocessing, we can guarantee that every feature has a constant number of 
neighboring features. This is how we can verify or update the closest feature pair in 
expected constant time. Fig. 5 shows how this algorithm works on a pair of simple 
convex polytopes. 

Sweep and Prune 

The exact polytope collision detection algorithm is fast enough to perform collision 
detection at real-time, in a multi-body environment of moderate size. However, the 
0(n 2 ) pairwise detection tests become a bottleneck for large n. By exploiting the 
spatial and temporal coherence, we use a method whose running time is a linear 
function of the number of objects and the number of pairs of objects in close vicinity 
of each other. 

For each object in the environment, we compute a smallest fitting, axis-aligned 
bounding box, which is dynamically resized using a hill-climbing algorithm based 
on the properties of convex polytopes and the principle of coherence. 

These bounding boxes are used to prune down the number of pairwise collision 
detection tests. We take the orthogonal projections of the bounding boxes onto the 
x, y, and z axes and sort the intervals in three separate lists. If the projections over- 
lap in all three dimensions, we activate the exact collision detection algorithm. At 
each frame, we sweep over the sorted lists of intervals, updating them and culling 
out the pairs of non-overlapping bounding boxes. 

Each polytope in the environment now has a bounding box. We only call the 
exact collision detection algorithm between two objects when their bounding boxes 
overlap. Due to coherence, the sorted lists of intervals, constructed from the projec- 
tions of the bounding boxes, change little from frame to frame, so we use a bubble 
or insertion sort to update these lists in expected linear time. Without the bound- 
ing box hierarchy, we would have to perform a quadratic number of pairwise exact 
detection tests. 
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Video Contents 



We have successfully demonstrated the algorithm in an architectural walkthrough 
environment and a multi-body simulation system. For simulations consisting of 
thousands of models, the overall algorithm has exhibited linear time performance 
in practice. 

We first demonstrate how the polytope collision detection algorithm [2] tracks 
the closest features between two convex polytopes in real time. We maintain and 
update the closest features using the Voronoi regions of the polytopes. Since the 
polytopes move only slightly between the frames, the temporal and geometric co- 
herence is well preserved. Thus, the expected running time of this algorithm is con- 
stant, independent of the complexity of the polytopes. 

This is followed by demonstrating the algorithm on an articulated body, namely 
a hand, represented by the union of multiple convex polytopes. To demonstrate the 
efficiency of the algorithm and effectiveness of the coherence, we keep track of all 
pairs of closest features between each finger and the nearby objects. However, using 
a subpart hierarchical tree representation as in [2] can eliminate the unnecessary 
computation. 

The over all collision detection algorithm is tested on a walkthrough and a multi- 
body dynamic simulation environment. The video footage of walkthrough was taken 
from a user's actual walkthrough experience. The collision detection algorithm 
performs very well in large environments composed of hundreds of polytopes. The 
kitchen, like most realistic environments, is fairly static; therefore, we can sort the 
bounding boxes of the polytopes in sub-linear time at each instance (this requires a 
modified sorting algorithm which sorts only the moving objects). We see very little 
degradation of frame rate when we add the capability of collision detection to our 
walkthrough environment. 

We further demonstrate a multi-body simulation with several complex objects 
in a moderately dense environment. The objects are placed in a bounded 3-dimensional 
simulation volume whose walls act as barriers. Unlike the walkthrough environ- 
ment, all the objects move independently. The dynamics of multi-body simulation 
assumes elastic collision. In order to show the asymptotic time bound on the algo- 
rithm, we increase the number of possible interactions quadratically, but the frame 
time using our collision detection algorithm increases only linearly. 
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The accompanying videotape illustrates two deterministic polynomial time methods 
for finding a separator for two nested convex polyhedra in 3-dimensions. 

1 Introduction 

Let Q C P be two convex polyhedra in R 3 . The problem of finding a separating 
convex polyhedron between P and Q with as few facets as possible is believed to be 
NP-haid [2], but one can find good approximations for it. Mitchell and Suri [4] cast 
the problem as a set cover problem. Let H = H(Q) denote the set of supporting 
hyperplanes of Q. Let the cover set system be (H, {H p : p £ dP}), where, for any 
point p, Hp consists of those hyperplanes in H for which p and Q lie on opposite 
sides. They show that a hitting set for the cover set system gives a subset of facets 
of Q whose supporting halfspaces define a polytope whose boundary lies between 
Q and P. 

A polytope Q' between Q and P such that each facet of Q' is supported by a 
facet of Q is called canonical. The smallest (with respect to the number of facets) 
canonical polytope contained in P can be obtained from a minimal hitting set of this 
set system and is within 3 times the optimal separating polyhedron (in the number 
of faces). Therefore, in this video, we just show how to obtain an optimal canonical 
separator. 



25 



7 



2 Video Contents 

We first familiarize the viewer with the two polyhedra Q C P; Q is seen inside 
by means of transparency (see Fig. 6, top). We explain how to obtain canonical 
separators, and how to construct the set system. The portion of the arrangement A 
of H(Q) on the boundary of P is visualized, and some rows of the set system are 
then shown sequentially, one corresponding to each facet of this arrangement. 

To explain the greedy method on this set system, an approach which is proposed 
by Mitchell and Suri [4], we show how the choice of a facet removes a portion of the 
boundary of P, and how the greedy method proceeds sequentially by choosing the 
facet which covers the biggest remaining portion (measured by the number of cells 
of A it separates from Q). After four iterations, we skip to the final greedy separator, 
which has 14 facets. In general, the greedy separator is guaranteed to have at most 
s op t log \ Q\ facets, where s opt is the size of an optimal canonical separator. 

Another algorithm is then animated. It is the weighting method described by 
Goodrich and the author in [1]. We show how putting a weight on the faces of Q 
induces a weight for the cells of the arrangement A, and how keeping only the heavy 
cells yields a separator which is nothing else than a (l/2c)-net for the planes H(Q) 
supporting faces of Q (see Fig. 6, bottom). We then explain the reweighting strategy 
in this context. Eventually, a separator is found that has 13 facets. In general, the 
weighted separator is guaranteed to have 0(s opt ) facets. 
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Figure 6: The top image shows two nested polyhedra. The 
bottom image shows the weighted algorithm in action. Specif- 
ically, the heavy cells are displayed in yellow on P, and the 
cover for those cells is a (l/2c)-netfor the planes with weight, 
displayed in green on the boundary of Q. 
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3 Technical Details 

Both algorithms are implemented and animated using GASP, a Geometric Anima- 
tion System created at Princeton University by Ayellet Tal and David Dobkin [3]. 
The implementation is written in C and runs on a Silicon Graphics Iris workstation. 
The algorithm is fully animated, in the sense that it produces an animation for any 
given data P and Q. The data used for the video were created with a random gener- 
ator of planes tangent to the sphere, then intersected to obtain Q, which was blown 
up to yield P. 
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The accompanying videotape shows an interactive visualization of weighted three- 
dimensional a-hullsfor static and dynamic spheres. The a-hull is analytically com- 
puted and represented by a triangulated mesh. The entire surface is computed and 
displayed in real-time at interactive rates. The weighted three-dimensional a-hulls 
are equivalent to smooth molecular surfaces of biochemistry. Biochemistry appli- 
cations of interactive computation and display of a-hulls or smooth molecular sur- 
faces are outlined. 

Introduction 

The a-hull has been defined as a generalization of the convex hull of point-sets by 
Edelsbrunner, Kirkpatrick, and Seidel [3, 4]. Given a set of points P, a ball b of 
radius a is defined as an empty a-ball if b n P = <f>. For 0 < a < oo, the a- 
hull of P is defined as the complement of the union of all empty a-balls [3]. It has 
been shown that it is possible to compute the a-hull of a set of points P from the 
Voronoi diagram of P. For a = oo the a-hull over the set of points P is the same 
as their convex hull. Edelsbrunner and Miicke [5] have defined the a-shape over 
P to be the polytope that approximates the a-hull over P by replacing circular arcs 
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of the a-hull by straight edges and spherical caps by triangles. An a-shape of a set 
of points P is a subset of the Delaunay triangulation of P. Edelsbrunner [3] has 
extended the concept of a-shapes to deal with weighted points (i.e. spheres with 
possibly unequal radii) in three dimensions. An a-shape of a set of weighted points 
P w is a subset of the regular triangulation of P w . 

The smooth molecular surface of a molecule is defined as the surface which an 
exterior probe-sphere touches as it is rolled over the spherical atoms of that molecule. 
This definition of a molecular surface was first proposed by Richards [7]. This sur- 
face is useful in studying the structure and interactions of proteins, in particular for 
attacking the protein-substrate docking problem. The analytic computation of the 
molecular surface was first done by Connolly [2]. 

Looking at the above definitions, one can see that the weighted a-hulls for three 
dimensions and the smooth molecular surface of a molecule with a probe-radius a 
have the same definitions. In this video we present the visualization of weighted 
a-hulls as used to define molecular surfaces. 

Our Approach 

Computation of a-hulls and a-shapes has traditionally been done by first construct- 
ing power diagrams or regular triangulations. Since these methods involve comput- 
ing the entire diagram or triangulation first and then culling away the parts that are 
not required, their complexity is 0(n 2 ) in time, where n is the number of points. 
This is worst-case optimal, since an a-shape in three dimensions could have a com- 
plexity of Q(n 2 ). 

However, when a-hulls are used as molecular surfaces, one can do better. Let 
M = {Si, . . . , S n }, be a set of spheres, where each sphere, Si, is expressed as a 
pair (ci,ri), Ci being the center of the sphere and n being the radius of the sphere. 
Collections of spheres representing molecules have two interesting properties: (i) 
the minimum distance dij between any two centers a and cj is greater than or equal 
to a positive constant l m in — the smallest bond-length in the molecule and (ii) the 
values of all the radii can be bounded from above and below by strictly positive 
values, 0 < r m i n < r« < r max . We take advantage of the first property to arrive 
at better running times for our algorithm. Stated simply, the first property says that 
the number of neighboring atoms within a fixed distance from any atom i, is always 
bounded from above by a constant k max that depends on the minimum spacing be- 
tween any two atoms. We refer to two atoms i and j as neighboring if it is possible 
to place a probe sphere such that it is contact with both Si and Sj. 
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Figure 7: The molecular surface ( alpha-hull) for crambin at 
probe-sphere radii (alpha value) of 1.4 Angstroms. 

A power cell (as described by Aurenhammer [1]) for a given sphere Si with 
respect to M can be computed as the intersection of n — 1 half spaces, each half- 
space contributed by the pair (Si, Sj), 1 < j < n, j ' ^ i. If the average number 
of neighbors for an atom is k, then for our purposes it is sufficient to just compute 
an approximation to the power cell, which we call a. feasible cell, as an intersection 
of k halfspaces (one halfspace contributed by each neighbor). This can be done 
in deterministic time 0(k log k). For n atoms, this task can be parallelized over n 
processors, each processor computing the feasible cell for one atom. To check if a 
feasible cell is non-null, we use a randomized linear programming algorithm that 
has linear expected time complexity and is quite fast in practice [8]. For details of 
our approach the interested reader can refer to [9] . 

We have used ideas from the theory of packing and covering of spheres to es- 
timate the value for k, the average number of neighboring atoms, for an atom in 
protein molecules. We can prove that for proteins, k < 140 for a probe-sphere ra- 
dius of 1.4A- the radius of a water molecule. Details about this result can be found 
in [10]. In practice we have found that for protein molecules and for a probe-radius 
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of 1.4 A, k is around 45. 

In the algorithms for computing the convex hull of a set of points, it is assumed 
that the points are in a general position, i.e. no more than d points lie on the same 
d—1 dimensional hyperplane. In reality this assumption often fails to hold, leading 
to problems. For example, planar benzene rings occur often in proteins, causing six 
carbon and six hydrogen atoms to be all coplanar. One of the recent approaches to 
solving this problem has been to perturb the input point set slightly to avoid these 
degeneracies. The approach of Emiris and Canny [6] perturbs the j th dimension of 
the i th point as: 

Pi,j( e ) = Pi,j + e {i 3 mod q)l < i < n, 1 < j < d (1) 

where e is a symbolic infinitesimal and q is the smallest prime greater than n. In- 
stead of performing exact integer arithmetic, we just perturb the centers of the spheres 
by the above scheme and that has worked quite well for us in practice. 

With no preprocessing, we can compute the molecular surface for a 396 atom 
Crambin and a probe-radius of 1.4 A, using 40 Intel i860 processors in 0.2 seconds. 
We have implemented a-hulls on Pixel-Planes 5, though the method is general enough 
to be easily implemented on any other parallel architecture. Our approach is analyt- 
ically exact. The only approximation is the degree of tessellation of spherical and 
toroidal patches. 

Applications in Biochemistry 

Interactive computation and display of molecular surfaces should benefit biochemists 
in three important ways. First, the ability to change the probe-radius interactively 
helps one study the surface. Second, it helps in visualizing the changing surface of a 
molecule as its atom positions are changed. These changes in atom positions could 
be due to user-defined forces as the user attempts to modify a molecular model on a 
computer. Third, it assists in incorporating the effects of the solvent into the overall 
potential energy computations during the interactive modifications of a molecule on 
a computer. 

Scope for further work 

At present we are not using any incremental temporal information in constructing 
these surfaces. Thus, if the atoms move slightly from their positions, the whole sur- 
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Figure 8: The molecular surface (alpha-hull) for crambin at 
probe-sphere radii (alpha value) of 10 Angstroms. 

face is recomputed from the beginning, although one could conceivably store the 
feasible cell description with each atom. Assuming the atoms of the molecule move 
along continuous trajectories, it should be possible to compute such surfaces (and 
indeed a-hulls and a-shapes) incrementally and efficiently by using the information 
from previous time steps. 
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